Numerical simulations of flow reversal in Rayleigh-Benard convection 
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We investigate numerically the statistical properties of the large scale flow in Rayleigh-Benard 
convection. By using an external random perturbation on the temperature field, we were able to 
decrease the effective Prandtl number of the flow while keeping the Rayleigh number relatively small; 
this increases the Reynolds number thus making possible the numerical investigation of the long- 
term flow statistics. We also propose a simple and quantitative explanation for the experimental 
findings on the statistical distribution of flow reversals and reorientations. 



Turbulent convection in a Rayleigh-Benard cell is char- 
acterized by the complex behaviour of the unstable ther- 
mal boundary layers and plumes which produce a large 
scale (turbulent) flow in the cell. Experimentally, it has 
been observed that the large scale "wind" may exhibit 
abrupt flow reversals [l[ and reorientations [2| > @] > > 
whose statistical properties have been the subject of some 
theoretical investigations [H , [H , Q . Up to now, there has 
not been any serious numerical simulations aimed at in- 
vestigating the statistical properties of large scale flow in 
Rayleigh-Benard systems. This is due to the fact that 
interesting statistical properties of large scale flow are ob- 
served (experimentally) only at relatively large Rayleigh 
numbers Ra, 10 8 — 10 11 , and after having collected data 
for rather long time periods (up to one year Q). From 
a numerical point of view, major computer facilities are 
needed in order to properly resolve large Ra conditions 
and simulating these flows for thousands of large eddy 
turnover times is, at present unfeasible. In this Letter, 
by using a suitable strategy, we show how to perform 
numerical simulations aimed at investigating the statis- 
tical properties of large scale flow in Raylaigh Benard 
convection while maintaining the computational require- 
ments at a reasonable level. The flow investigated in 
this Letter is that developing in a cylindrical cell of as- 
pect ratio (diameter d over cell height h) Y = d/h = 1 
heated from below and cooled from above with an adi- 
abatic side wall. All the surfaces boundary conditions 
are no-slip. 10 azimuthally equi-spaced numerical ideal 
probes are located on a circle halfway between the plates 



at a radius r p = 0.2h; the "probes" provide simultaneous 
point-wise measurements of temperature and of the three 
velocity and vorticity components. The flow is solved by 
numerically integrating the three dimensional unsteady 
Navier-Stokes equations with the Boussinesq approxima- 
tion: 
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where D/ Dt = dt + u • V , z the unity vector pointing 
in the opposite direction with respect to gravity, u the 
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FIG. 1: Time evolution of the vertical velocity in one ideal 
probe. Many rather abrupt changes in the signal are observed. 
The numerical simulation are performed using (I I 121) with an 
external random forcing on the temperature equation ([2} . In 
the inset we show the same quantity for the same value of 
Ra = 6 x 10 5 without the random perturbation on the tem- 
perature field. 



velocity vector, p the pressure (separated from its hydro- 
static contribution) and 9 the non dimensional temper- 
ature. The equations have been made non dimensional 
using the free-fall velocity U = y/gaSE, the distance 
between hot and cold plates h and their temperature dif- 
ference A = T/j — T c ; the non-dimensional temperature 9 
is defined 9 = (T - T c )/A so that < 9 < 1. The above 
equations have been written in a cylindrical coordinate 
frame and discretized on a staggered mesh by central 
second-order accurate finite-difference approximations; 
the resulting discretized system is solved by a fractional- 
step procedure with the elliptic equation inverted using 
trigonometric expansions in the azimuthal direction and 
the FISHPACK package [| for the other two directions. 
The numerical method is the same as that described in 
fl0| and [ll| where further details of the numerical pro- 
cedure can be found. The numerical experiments were 
performed at Pr = 0.7 and Ra = 6 x 10 5 on a grid of 
33 x 49 x 97 nodes, respectively, in the azimuthal, radial 
and vertical directions; this grid was found to be suffi- 
cient in pj| for the used flow parameters. Owing to the 
relatively low value of Ra, large scale flow reversals or 
reorentations are expected to occur on extremely long 
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FIG. 2: Time evolution of the angular orientation O of the 
interface between warm rising fluid and cold falling fluid. In 
the inset we show the probability distribution of the vertical 
velocity as computed in one single probe. 



time scales, too long to achieve an acceptable statistical 
convergence within a reasonable computational effort. In 
order to increase the turbulent fluctuations, one can ei- 
ther increase the value of Ra or alternatively decrease 
the Prandtl number (which is equivalent to increase the 
Reynolds Re number of the system), both methods im- 
plying an increase of numerical resolution. Alternatively, 
we can artificially increase thermal fluctuations by adding 
to the right hand side of an extra forcing term /, 
6 — correlated in space and time. More precisely the 
explicit expression is / = evA?0 where e is the ampli- 
tude of the perturbation, At is the time step size and (f> 
is a space dependent white-noise random number. The 
physical motivations for introducing / are twofold: first 
of all, it has been pointed out (Q) that small scale turbu- 
lent fluctuations can be quantitatively considered as an 
external random forcing acting on the large scale flow. 
Thus, one may consider / as the overall effect of small 
"unresolved" turbulent motion for a smaller value of Pr. 
Secondly, one can think that the effect of the external 
noise is to increase thermal diffusivity and, as a conse- 
quence, to increase Re at fixed Nussclt number Nu Both 
arguments, however, have the same physical meaning. 
In this particular example the Nusselt number resulted 
Nu = 4.1 ± 0.1 with and without the external noise on 
the temperature field, as expected, because the Nusselt is 
essentially independent of Pr, i.e. on thermal diffusivity. 
In figure {1} we show the vertical velocity at one single 
ideal location as function of time for e = 0.1(T/j — T c ), 
while in the inset we show the same result for e = 0. As 
we can see, a rather clear increase of the number of re- 
versals is observed by using the external random forcing 
/. The value of e ~ 0.1 (Th — T c ) has been chosen by trial 
and error, i.e. not too strong with respect to the deter- 
ministic dynamics and not too small to be irrelevant for 
the time scale of the reversal. As a consistency check, 
however, we have run similar simulations using half and 
twice the above value of e obtaining similar statistics for 
the large scale flow. The main simulation was run for 



10 5 time units that on account of the cell geometry cor- 
respond to about 10 5 /-7r w 3 x 10 4 large eddy turnover 
times. Graphical inspection shows that the flow in the 
cell is roughly divided into two halves, one with a rising 
warm current (positive vertical velocity) and the other 
with cold sinking fluid (negative vertical velocity). The 
two regions are separated by an "interface" where the 
vertical velocity is close to zero. Because of the cylindri- 
cal geometry, the interface has no preferred orientation 
and, therefore, it is fluctuates in time. By using the same 
procedure described in Q and 0] we can compute the 
instantaneous orientation of the interface. As a con- 
sistency check we have applied the same procedure using 
either the temperature and the vertical velocity signals of 
the probes, always obtaining the same results. In figure 
(|2"|). we show the time behaviour of 0, while in the inset 
we show the probability distribution of the vertical ve- 
locity as obtained by measuring it on a single probe. The 
behaviour of is in close agreement with the experimen- 
tal data of [i[ and [H, while the probability distribution 
of the vertical velocity shows a bimodal shape in agree- 
ment with [l|. Next we can study the statistical proper- 
ties of our large scale flow. In particular, the existence 
of a bimodal distribution for the vertical velocity on a 
single probe, allows us to study the statistical proper- 
ties of the (random) switching time (r) between the two 
states identified as the maxima of the probability shown 
in ([2]). More precisely, we define r as the time interval 
between two successive zero-crossings of the vertical ve- 
locity at one ideal probe and we compute its probabity 
distribution P(r). This study was first performed by [l[ 
where an estimate of the vertical velocity was obtained 
by using the correlation time of two close temperature 
probes. In [l[ and P(t) was found to be a power 
low for small enough r, while at large r an exponential 
cutoff in P(t) was observed. This is a rather peculiar 
result which is difficult to explain by using simple ar- 
guments based on the theory of stochastic differential 
equations or chaotic dynamics Q. Our numerical simu- 
lations agree with the experimental results as shown in 
figure ([31 circles). Following Q we also computed the 
probability distribution Pg of Tg defined as the time be- 
tween two events where both <50 = Q(t + 5t) — 0(i) and 
dQ/dt = 50/ 5t arc larger than given thresholds; we have 
used for the thresholds the same values as [1] and have 
verified that the statistical results are robust with re- 
spect to the arbitrarily selected values. In the inset of 
figure [3] we show the log-lin plot of Pg and compare it 
with an exponential distribution (line). As observed in 
the experimental analysis of [|[ , the probability Pg seems 
to be fitted rather accurately by an exponential distribu- 
tion. Figures and ([3]) show our main results, i.e. our 
numerical simulations reproduce qualitatively and quan- 
titatively the experimental findings. Next, we want to 
discuss briefly some theoretical ideas which may help in 
explaining the power law behaviour of the probability 
distribution P(t) as computed by [l[ and reproduced in 
our simulations. From a theoretical point of view, we 
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FIG. 3: Circles:log log plot of the probability distribution 
P(t) of the random switching time r. Note that P shows a 
power law behaviour shown by a straight line in the figure. 
Inset: log-linplot of the probability distribution Ps(ts) of 
Tg following [31 . Note that the best fit is consistent with an 
exponential shape. 
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FIG. 4: Probability distribution of log(r) for the ipo defined 
in Q (line) as compared with the experimental data (circles) 
and the case of an overdamped brownian particle in a double 
well potential (see text for a description). 




FIG. 5: Log-log plot of the probability distribution P(r) of 
the random switching times r between the maxima of w. 
Small black circles are obtained using ([6jl with O a random 
noise correlated in time. Line with symbol refers to the exper- 
imental data by In the insert, we show the same quantity 
computed with O obtained by our numerical simulations, see 
figure ©. 



can state the problem as follows: Let Q be a random 
variable whose probability distribution is bimodal with 
maxima in ±Qo an d let tq be the (random) switching 
time between the maxima. Is there any "simple" model 
for which P(tq) shows a power law behaviour for small 
tq? We will show that such a model does exist and it 
may be relevant for thermal convection. To this aim, let 
us consider the Landau- Ginzburg equation (LGE) in two 
space dimensions for the field ip(x, y, t): 



d t tjj = mip — gijj 3 + vA<ip + crr)(x, y, t) 



(3) 



where to, <?, v, a are real positive variables and 77 is a white 
noise delta correlated in space and time. We assume ip 
defined on a periodic box of finite size L and we indi- 
cate with (..) s = L~ 2 J ..dxdy the space average. We can 
decompose ip in its space average ipo = (ip) s and a fluc- 
tuation 0, i.e. ip = ipo + (f>. Then the equation for ipo 
reads: 

d t ipo = (m - 3.g((/) 2 ) s )Vo - <#o + ( 4 ) 

where rj is a white noise in time. More precisely, let 
us consider ([3]) defined on a regular lattice of N x N 
points of spacing Ax = Ay = L/N. Then a s = a/N 2 
and we assume a s = const, for any N (see [HI for a 
correct definition of the problem). It is worth remind- 
ing that the equilibrium probability distribution of ([3]) 
has been widely used as a model of second order phase 
transition in the infinite volume limit. For a suitable 
choice of the parameters m,g and v, the term (4> 2 ) s can be 
large when ipo becomes small, i.e. the system at ipo ~ 
shows strong fluctuations as it is usual in second order 
phase transitions. In this case, the linear term in (|4|). 
i.e. (to — 3g(<p 2 ) s ), can occasionally change sign and the 
switching between the two maxima of the probability dis- 
tribution of ipo can be more frequent. This is equivalent 
to say that the switching time is no longer controlled 
by the mechanism leading to an exponential distribution 
and one can observe a scale (in time) invariant proba- 
bility distribution, i.e. a power law. In figure (j4|) we 
show the probability distribution of the switching time 
of ipo obtained by numerically integrating equation l[5]). 
for N = 32, v = 0.1, to = 1, g = 9 and a s = 0.05. In 
particular we plot the probability distribution of log(r), 
where t is the random switching time computed by us- 
ing the parameter ipo- Note that if P(t) ~ r _1 then 
P(log(r)) ~ const. In the same figure we show the prob- 
ability distribution of log{r) for the experimental data by 
[H) and for the well known case of an over-damped Brow- 
nian particle in a double well potential with the same pa- 
rameter to and g and stochastically perturbed by a white 
noise with variance a s . In the last case, we know that the 
probability distribution of r is exponential and, indeed, 
in figure (|1]), P{log(r)) for the double well potential is 
not constant. On the other hand, P{log(T)) is constant 
for almost one decade for the numerical simulation of ([3]) 
using ipo as order parameter. Thus, we may tentatively 
state that the equation and in particular ([?]), can be 
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considered as a conceptual model for the experimental 
results of From a physical point of view, we remark 
that there exists a non trivial feedback between the vari- 
ance of the fluctuations in space and the value of the 
order parameter ipQ. It is now tempting to understand 
whether the same "physical" mechanism can be observed 
in our numerical simulation of the Raylcigh-Benard tur- 
bulence. According to our previous discussion, the over- 
all picture which emerges is quite simple and agrees with 
the original description of [J]. As a first approximation 
we can think that the large scale flow rotates randomly 
in time with some orientation Q(t) performing a random 
walk in the interval [0, 2tt]. The vertical velocity in a 
single point is almost insensitive to unless the inter- 
face is close to the ideal probe. A suitable model for 
0(i) can be build as follows. Let us consider two vari- 
ables x\ = r cos((9) and x 2 = 7'sin(#) which satisfy the 
following equation 

dxi dV . , . 

— = h noise (5) 

at oxi 

where V = — l/2r 2 {R 2 — r 2 ) and r 2 = x\ + x 2 , i.e. r ~ R 
for most of the time. We can imagine that 9 is precisely 
the variable describing the interface direction of the 
large scale wind. Then denoting by w the value of the 
vertical velocity, a qualitative model of w as a function 
of 6{t) is simply given by: 

x 2 sin(6) 

W = ; = ; 6 

y/x 2 + A y/sin 2 {6)+A/r 2 

where the parameter A is a measure of the interface thick- 
ness, whose fluctuations are controlled by r. Using © 
and our simple model, one can easily show that: 

dW n 

— = aw — (a + bA)w + noise (7) 
dt 



where a = R 2 — 2x\ and b = 2.. Equation ([7]) has a 
form of an over-damped Brownian particle in a double 
well potential. However the term a can change sign (ex- 
actly as in @) because of the (strong) fluctuations of X\. 
In order to compute the probability distribution of the 
switching time r between the two maxima in the proba- 
bility of w, we have numerically integrated equation (|5j). 
with R — 1, A = 0.1 and the variance of the noise equal 
to 0.1. The final result is shown in figure §5§ (small black 
circles) and compared with the experimental finding of 
[l[ (line with symbol). As we can sec P{t) is close to 
a power law with slope —1 in agreement with the ex- 
perimental findings of [![. In the same figure, we show 
the probability distribution P(t) computed by using the 
approximate value w ~ sin (&)/ v / sin 2 (e) + A and the 
value of obtained in the numerical simulation (see fig- 
ure (J2J)). A good agreement is observed with a power 
law P(t) ~ t _1 . Thus © seems to capture the physi- 
cal reason for the observed numerical and experimental 
results in the statistical properties of flow reversal. In 
summary, we have shown that the statistical properties 
of large scale flow in Rayleigh Bcnard convection can be 
efficiently simulated by using an external random forc- 
ing in the temperature field. The effect of the forcing 
is to increase thermal diffusivity, i.e. to decrease the 
Prandtl number. Consequently, the effective Reynolds 
number achieved in this way is larger with respect to 
the one without external forcing. Our numerical simu- 
lations agrees qualitatively and quantitatively with the 
experimental findings of [l[ , [|[ and 0] . Finally, we have 
discussed some theoretical ideas which may be able to 
explain the statistical properties of flow reversals. 
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